Using the Wigner-Ibach Surmise to Analyze Terrace- Width Distributions: History, 

User's Guide, and Advances 
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A history is given of the applications of the simple expression generalized from the surmise by 
Wigner and also by Ibach to extract the strength of the interaction between steps on a vicinal 
surface, via the terrace width distribution (TWD). A concise guide for use with experiments and a 
summary of some recent extensions are provided. 

PACS numbers: 68.35.Md, 05.40.-a, 68.37.Ef 



Dedicated to Prof. Harold Ibach, with profound grati- 
tude, on the occasion of his retirement 

INTRODUCTION 

Measurement of the terrace width distribution (TWD) 
P(£) of vicinal surfaces is now used routinely to find 
the dimensionless strength A of the elastic repulsion be- 
tween steps. The University of Maryland [l| and the 
Forschungszentrum-Jiilich Q have been at the vanguard 
of this progress. Use of an extension of the Wigner sur- 
mise from random matrix theory has resolved ambiguities 
on how best to estimate A from the variance. This paper 
discusses the history of this development and the crucial 
role played by Harald Ibach in divining from physical 
insight an analytic expression for the TWD associated 
with the special case A = that is essentially the same 
as the simple expression that Wigner surmised describes 
(albeit, ultimately, not exactly) the distribution of the 
energy differences of adjacent levels in a nucleus when the 
coupling has unitary symmetry. This insight spawned my 
long collaborative effort with FZ-Jiilich to develop gen- 
eralizations appropriate to arbitrary vicinal surfaces and 
to corroborate the viability of the resulting formulation 
to account for extensive experimental data. In addition 
to a personal history of this progress (with a few new re- 
sults) , I collect (with enhancements) in one place several 
tables from publications and present a concise "User's 
Guide" for applying the formalism and a short discus- 
sion of some recent developments, e.g., for azimuthally 
misoriented vicinal surfaces and for non-equilibrium sit- 
uations. Space limitations preclude fuller discussions or 
an authoritative update on experimental progress. Three 
of the four figures are heretofore unpublished. Also in- 
terwoven are comments stimulated by several penetrat- 
ing questions, many by Harald Ibach, that arose during 
talks I gave in Germany that forced me to sharpen my 
thinking. These are topics typically skirted in publica- 
tions as being well-established, but my experience is that 
they are not generally well understood. The use of "we" 
below is not a polite affectation but instead reflects the 
crucial role played by various of my many collaborators 
in this extended series of studies. 



FIG. 1: Plots of generic lattice configurations and associated 
TWDs, to illustrate essential features, as discussed in text, 
of steps on a vicinal surface (with no energetic interactions) . 
Top left: "perfect" cleaved crystal. Top right: straight steps 
placed randomly, corresponding to a strictly ID model. Lower 
left: meandering steps in a TSK model. Lower right: TWDs 
P(s) associated with these distributions. The downward ar- 
row emphasizes the greatly decreased chance of finding steps 
at very small separation when they meander, as compared to 
straight steps, due to the entropic repulsion. 



To set the stage (cf. Fig. [T]) , the direction perpendic- 
ular to the terraces (which are densely-packed facets) is 
typically called z. In "Maryland notation" the normal to 
the vicinal surface lies in the x — z plane, and the dis- 
tance £ between steps is measured along x , while the steps 
run along the y direction. In the simplest and usual ap- 
proximation, the repulsions between adjacent steps arise 
from two sources: an entropic or steric interaction due to 
the physical condition that the steps cannot cross, since 
overhangs cannot occur in nature. The second comes 
from elastic dipole moments due to local atomic relax- 
ation around each step, leading to frustrated lateral re- 
laxation of atoms on the terrace plane between two steps. 
Both interactions are proportional to l/£ 2 . 

To illustrate the essence of P(£), we consider in Fig. [1] 
its shape for 3 idealized configurations. There is just one 
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characteristic length in the x direction, namely the av- 
erage step separation (£). (Contrary to widespread mis- 
conception, (I) need not be a multiple of, or even sim- 
ply related to, the substrate lattice spacing. It is the 
step height times cot 0, where </> is the arbitrary mis- 
orientation, as in shown in Fig. [1]) Therefore, we con- 
sider P(s) = (£)~ 1 P(£), where s = 1/(1), a dimensionless 
length. For a "perfect" cleaved crystal, P(s) is just a 
spike S(s — 1). If, as intrinsic to ID models, the steps 
are imagined as uncooked spaghetti dropped at any po- 
sition with probability !/(£), P(s) is a Poisson distribu- 
tion exp(— s). Actual steps do meander, as one study 
most simply in a terrace-step-kink (TSK) model. In this 
model, the only excitations are kinks (with energy e) 
along the step. (This is a good approximation at low 
temperature T since adatoms or vacancies on the terrace 
cost several e [4e in the case of a simple cubic lattice]. 
The entropic repulsion due to step meandering dramat- 
ically decreases the probability of finding adjacent steps 
at £ <C (£)■ To preserve the mean of one, P(s) must also 
be smaller than exp(— s) for large s. 

If there is an additional energetic repulsion A/£ 2 , the 
magnitude of the step meandering will decrease, narrow- 
ing P(s). As A — > oo, the width approaches a delta func- 
tion. Note that the energetic and entropic interactions 
do not simply add. In particular, there is no negative (at- 
tractive) value of A at which the two cancel each other. 
Thus, for strong repulsions, steps rarely come close, so 
the entropic interaction plays a smaller role, while for 
A < 0, the entropic contribution increases, as illustrated 
in Fig. [2] and worked out in detail below. 



HISTORY 

Investigation of the interaction between steps on vic- 
inal surfaces is a core part of the flourishing field of ex- 
ploring the properties of these technologically important 
and scientifically rich systems, as discussed in several ex- 
cellent reviews The earliest studies seeking to 
extract A from TWDs used the mean-field-like Gruber- 
Mullins Q approximation, in which a single active step 
fluctuates between two fixed straight steps 2(£) apart. 
Then the energy associating with the fluctuations x(y, t) 
is 



AS = -/9(0)L„ 



mv))\n 



dx 

dy 



dy, (l) 



where /? is the step free energy per length (or line tension 
J5| ) for a step at orientation 9 relative to the mean direc- 
tion of the step (and the direction of the fixed, bounding 
steps), and L y is the size of the system along the mean 
step direction (i.e. the step length with no kinks). We ex- 
pand 13(6) as the Taylor series /3(0) + j3'(0)6 + y 2 /3"(0)6> 2 
and recognize that the length of the line segment has 




FIG. 2: Illustration of how entropic repulsion and energetic 
interactions combine, plotted vs. the dimensionless energetic 
interaction strength A = a/?/(feT) 2 . The dashed straight 
line is just A. The solid curve above it is the combined en- 
tropic and energetic interactions, labeled A e s for reasons ex- 
plained below. The difference between the two curves at any 
value of the abscissa is the dimensionless entropic repulsion for 
that A. The decreasing curve, scaled on the right ordinate, is 
the ratio of this entropic repulsion to the total dimensionless 
repulsion A e s- It falls monotonically with A, passing through 
unity at A — 0. See the discussion accompanying Eq. (|16p for 
more information and explicit expressions for the curves. 



increased from dy to dy/ cos 9 w dy(\ + 1 /2 9 2 ). For close- 
packed steps, for which /3'(0) = 0, it is well known that 
(using 9 w tan 9 = dx/dy), 



AS 



/3(0) 



dx 
dy 



dy, /3(0) = /?(0)+/3"(0), (2) 



where j3 is the step stiffness [6J. Most treatments gloss 
over the fact that the stiffness has the same definition for 
steps with arbitrary in-plane orientation. The key point 
is that to create such steps, one must apply a "torque" 
which exactly cancels the linear term (3' (0)9. Stasevich 
8] provides a more formal proof. 

Since x(y) is taken to be a single- valued function that is 
defined over the whole domain of y, the 2D configuration 
of the step can be viewed as the worldline of a particle 
in ID by recognizing y as a time-like variable. Since the 
steps cannot cross, these particles can be described as ID 
fermions. We also see from Eq. J2]) that in the (1+1)D 
formulation, dx/dy — ► dx/dt, with the form of a velocity, 
so that the stiffness plays the role of a mass; indeed, it 
serves as the inertial parameter of steps in this fermion 
perspective (though not with regard to actual dynamics 
in response to external forces Q). Moreover, stiffness is 
one of the three ingredients of the very-successful step- 
continuum model 

Pursuing this analogy for polymers in 2D, de Gcnncs 
[Io| showed nearly 4 decades ago that this problem could 
be mapped into the Schrodinger equation in ID, with the 
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thermal energy k B T replacing h. Then in the Gruber- 
Mullins approximation 0], the step with no energetic 
interactions becomes a particle in a ID infinite well of 
width 2{£), with well-known groundstate properties 



Me) * s[n {w)) 5 p{s) = sin2 ( t) 5 E ° 



{-Kk B Tf 



8/3 (£) 2 
(3) 

Thus, it is the kinetic energy of the ground state in the 
quantum model that corresponds to the entropic repul- 
sion (per length) of the step. In the exact solution for 
the free energy expansion of the equilibrium crystal shape 
[T]| , the numerical coefficient in the corresponding term 
is 1 /e rather than Note that P(s) peaks at s = 1 and 
vanishes for s > 2. 

Suppose, next, that there is an energetic repulsion 
U{£) = A/£ 2 between steps. In the ID Schrodinger equa- 
tion, the prefactor of -d 2 ip(t)/d£ 2 is (fc B T) 2 /2/3; hence, 
A only enters the problem in the dimcnsionless combi- 
nation A = Af3/(k B T) 2 [12]. In the Gruber-Mullins pic- 
ture, the potential (per length) experienced by the single 
active particle is (with I = I — {£)): 



U(£) = 



A 



A 



(£-(£)) 2 {£+{£)? 



2A 

W 



6A£ 2 10A£ 4 



(4) 

The first term is just a constant shift in the energy. For 
A sufficiently large, the particle is confined to a region 
|^| <C (£), so that we can neglect the fixed walls and the 
quartic term, reducing the problem to the familiar simple 
harmonic oscillator, with the solution: 



^oW«e- ? / fa! ; P G (s) 



1 



(TgV 2n 



exp 



2a 2 G 



(5) 

where aa = (48A)~ 1/4 and w = a G (£). 

For A of or 2, it turns out, as we shall see, that the 
TWD can be computed exactly. For these cases, Eqs. (O 
and |[SJ), respectively, provide serviceable approximations. 
It is Eq. ([5]) that is prescribed for analyzing TWDs in 
the most-cited resource on vicinal surfaces [l|. Indeed, 
it formed the basis of initial successful analyses of exper- 
imental STM data [13j. However, it has some notable 
shortcomings. Perhaps most obviously, it is useless for 
small but not vanishing A, for which the TWD is highly 
skewed, not resembling a Gaussian, and the peak, corre- 
spondingly, is significantly below the mean spacing. For 
large values of A, it significantly underestimates the vari- 
ance or, equivalently, the value of A one extracts from the 
experimental TWD width: Ihle et al. 14j point out that 



in the Gruber-Mullins approximation the TWD variance 
is the same as that of the active step, since the neigh- 
boring step is straight. For large A the fluctuations of 
the individual steps on an actual vicinal surface become 
relatively independent, so the variance of the TWD is 



the sum of the variance of each, i.e. twice the step vari- 
ance. Given the great (quartic) sensitivity of A to the 
TWD width, this is problematic. As experimentalists 
acquired more high-quality TWD data, other approxi- 
mation schemes were proposed, all producing Gaussian 
distributions with widths cx yl -1 / 4 , but with proportion- 
ality constants notably larger than 48~ 4 / 4 = 0.38. 

For the "free-fermion" (A = 0) case, Joos et al. [TBI ] 
developed at the beginning of the 1990's a sequence of 
analytic app roximants to the exact but formidable ex- 



pression [16|, |l_Jj for P(s). The procedure is based on a 
discrete version of the problem, representing the fermion- 
like steps in terms of second-quantized operators and 
taking note that the TWD is not just a pair-wise step- 
step correlation function but actually a many-particle 
correlation function in which we demand that no step 
lie between the steps at and £. The nth approxi- 
mant behaves well at smaller s but eventually passes 
through around s = n and then behaved non-physically. 
Specifically, the leading behavior of each approximant is 
1 — [sin(7rs)/(7rs)] 2 i=s Y3(7rs) 2 . In the asymptotic region, 
a different approach shows leading behavior for large s 
is P(s) cx s 7 / 4 exp[— y$(ns) 2 ], where the proportionality 
constant has the numerical value ^9.46. (The analytic 
form is in Ref. EH|.) This expression works to surpris- 
ingly small s, reproducing the peak semi-quantitatively. 
(It would have gone to even smaller s than indicated in 
Fig. 4 of Ref. [15| had the two leading as ymp totic correc- 
tions not been included!) Finally, Ref. [15| as well as a 
slightly earlier paper [l8( draw the analogy between the 
TWD of vicinal surfaces and the distribution of spac- 
ings between interacting (spinless) fermions on a ring, 
the Calogero-Sutherland model [3, [20| . which in turn 
for three particular values of the interaction — in one case 
repulsive (A = 2), in another attractive (A = — Yt), and 
lastly the free-fermion case (A — 0) — could, astonish- 
ingly, be solved exactly by connecting to random matrix 
theory [l7|, UK; Fi S- 5 01 R ei - depicts the three re- 
sulting TWDs. 

This was the state of affairs when Harald Ibach pre- 
sented me (seemingly sometime during his sabbatical stay 
in College Park during 1996-7) with Fig. [3l in which he 
achieved outstanding agreement with [numerical approxi- 
mation visually indistinguishable from] the exact solution 
by plotting what essentially amounts to 



Pi(s) = ^|s 2 exp(--s 2 ), 



(6) 



where the subscript of P refers to the exponent of s, 
and the pair of numerical factors ensure that the distri- 
bution is normalized and has a mean of one. He chal- 
lenged me to explain his insight, and Fig. [3] graced my 
office wall, greeting me each morning, for well over a 
year. Not hearing back from me, he put the remarkable 
expression (albeit with some typos that obscured its po- 
tency) into Ref. On Feb. 2, 1998, Safi Bahcall gave 
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FIG. 3: The graph for the expression deduced by H. Ibach for 
the TWD of steps with no energetic repulsion (solid curve), 
essentially Eq. surmised by Wigner, while the points are 
generated using the formalism for the approximants, given in 
Ref. El 



an intriguing condensed-matter seminar at University of 
Maryland on "Superconductivity and Random Matrices" 
(cf. Ref. [1^]). In subsequent discussions he rifled through 
a preprint of Guhr et a/.'s seminal review of random ma- 
trix theory [24]], and I spotted Fig. 12 therein, noting the 
similarity to Fig. 5 in Ref. as well as the accompa- 
nying discussion of the Wigner surmise, their Eq. (3.50) 
being the first I had seen of what will be Eq. (|9|) below. 

The other two curves correspond to corresponding 
Wigner surmises for those cases: 



/ 6 4 2 
s exp -— s 
9ir 



p i( s ) = 2 sl(5x P 



(7) 



(8) 



In random matrix literature, the exponent of s, viz. 1, 
2, or 4, is called /?, due to an analogy with inverse tem- 
perature in one justification. To avoid possible confusion 
with the step free energy per length {3 or the stiffness 
[3 for vicinal surfaces, we have called it instead by the 
Greek symbol that looked most similar, g. The Appendix 
offers transparent arguments on how the three kinds of 
symmetry lead to the associate exponents 1, 2, and 4. 
As seen most clearly by explicit plots, e.g. Fig. 4.2a of 
Haake's text [25[, Pi(s), P%(s), and P^s) are excellent 
approximations of the exact results for orthogonal, uni- 
tary, and symplectic ensembles, respectively, and these 
simple expressions are routinely used when confronting 



experimental data in a broad range of physical problems 
[24 HBJ . (The agreement is particularly outstanding for 
P2(s) and i-4(s), which are the germane cases for vicinal 
surfaces: the variance is 1% below and 0.4% above the 
exact values, respectively; for Pi(s) it is 4-1/2% below. 
This agreement is significantly better than any other ap- 
proximation (cf. Table 2 of Ref. [2||) and far better than 
the Gruber-Mullins approximation, as depicted in Fig. 1 
of Ref. [13, Fig. 1 of Ref. [H, and Fig. 2 of Ref. 

With the values g — 1, 2, and 4. the three specific 
expressions Eqs . (JHHE1) comprising the so-called Wigner 



surmise 



24|, [25j can be written as a single formula 



P e (s) = a e s e exp I 



-b e s 2 



(9) 



where the constants b g , which fixes its mean at unity, and 
a e , which normalizes P{s), are 



' g+2 \ 
, 2 1 



, 2 / 



2[F(^) ] e+1 _26 



<>+l)/2 



e+2 



(10) 

As noted above the Calogero-Sutherland model pro- 
vides a connection between random matrix results, no- 
tably the Wigner surmise, and the distribution of spac- 
ings between fermions in ID interacting with dimension- 
less strength A. Specifically, 



.4 



(I- 1 ) 



= 1 + V1 + 4A 



(11) 



For an arbitrary system, there is no reason that A should 
take on one of the three special values. However, it seems 
impossible to generalize the arguments in the Appendix 
to general values. For our purposes, the obvious solu- 
tion is to use Eq. (jlip for arbitrary g or A. Curiously, 
this form has not been applied in conventional investi- 
gations involving RMT, but instead other analytic phe- 
nomenological expressions, e.g. those proposed by Brody 
0] and by Izrailev [H are used (cf. Ref. Q). Such ap- 
plications typically involves mixtures of systems of two of 
the symmetries; neither they nor other models involving 
crossover between well-defined symmetries 36] are analo- 



gous to systems based on the Calogero-Sutherland model 
for arbitrary g. 

For arbitrary g there is no symmetry-based justifica- 
tion of distribution based on the Wigner surmise of Eq. 
©. Nonetheless, we have argued that it provides a 
viable, arguably optimal interpolation scheme between 
the two special values of g and also out to the Greno- 
ble expression for nearly-infinite repulsion; 27, [28j we 
have also used it successfully to analyze experimental 
data. [28, 37, 38| For brevity, we refer hereafter to this set 
of formulas, Eqs. (|9I10I) as the GWD (generalized Wigner 
surmise); elsewhere |26l. I28I [3e| we have called it CGWD 
(continuum generalized Wigner distribution). 

The moments of the GWD can be expressed simply in 
terms of b n : 
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TABLE I: Tabulation of various measurable properties of terrace-width distributions P(s) [where s is the terrace width nor- 
malized by its average value] based on exact results at the three soluble values of the dimensionless interaction strength A, the 
corresponding generalized Wigner distribution (GWD) expression, and the various Gaussian approximations: Gruber-Mullins 
(GM), modified Grenoble [3, H3, S3 (G), and Saclay IH(S). (In the original Grenoble approximation 0, a 2 tx A~ 1/2 
rather than A eS = 2/g, but with the same prefactor as indicated in this table.) SHO =>■ simple harmonic oscillator, i.e. 
uniformly spaced steps [energies] . The extreme case g=Q, for which exact results are trivial, is included to dramatize trends in g 
[3^ |. As g increases, the TWD becomes narrower, more symmetric, and more nearly Gaussian. Anticorrelations of neighboring 
terrace-width fluctuations increase. For the three exactly-solvable [non-trivial] cases, the GWD provides an excellent approx- 
imation, far better than any alternative. [Adapted from Ref. [2^], with most results for the "Exact" case from Refs. fl7l |24|. 
with new evaluations for symplectic case from Ref. [33j|.] 



Property Case 


g — g — 1 £) = 2 {3 = 4 £>^oo 
"Random" Attractive Non-interact Repulsive Extreme rpl. 
4 = 0" A = —1/4 1 = i = 2 i = £ 2 /4 


i=f(f-l) A eS = (|) 2 =i+| 


Underlying TL symmetry 


[Poisson] orthogonal unitary symplectic [SHO+phonons] 


Variance Exact 
a 2 = fi 2 GWD: [(g + 1) /2b e ] - 1 
= /4 - 1 GM (all): 05/87^ = 0.139i" 1/2 
GM (NN): 1/V4S = 0.144i" 1/2 


1 0.286 0.180 0.1041 0+ 
0.5708 0.2732 0.1781 0.1045 O^OOfT 1 
0.1307 0.0981 0.278e _1 
0.1307 0.1021 O.2890" 1 


Alternative G (all): 1 f' 27I }Z cos ^ =0.247 A~ 1/2 

V ' TV Jo <p(27r — 4>) cn 

estimate G (NN): y / 2/3n 2 =0.260A~^ 2 
of a 2 S: 2/tt 2 = 0.203ij ff 1/2 


0.1747 0.495£ -1 
0.1838 0.5202- 1 
0.203 0.101 O.4O50" 1 


Neighboring Exact cov a (si,S2) 
terraces Exact ((si + s 2 — 2) 2 ) 


-0.27 -0.31 -0.34 

2 0.416 0.248 0.138 0+ 


Peak Exact 
position GWD: (g/2b e ) 1/2 
Smax Gruber-Mullins 


0.77 0.8840 0.941 1" 
0.7979 0.8862 0.9400 1 - O^SOg' 1 
1 1 1 


Skewness Exact 

= GWD: /4 = (q + 2)/2b B 
(fj.3 — l)/cr 3 — 3/cr Gaussian Approx. 


2 0.4972 0.350(1) 
0.9953 0.6311 0.4857 0.3542 0.707£~ 1/2 




Kurtosis Exact 
M4/cr 4 GWD: fi'Ja 4 = {g + 3)/(g+l) 
Gruber-Mullins 


9 3.1 3.027(1) 
3.8691 3.2451 3.1082 3.0370 3 + 0.750£T 2 
2.4062 3 3 



"Covariance of adjacent terrace widths: cov(si, s 2 ) = -. <(si ~ (si>)(s2 ~ (i,2>)> 1l/2 = a- 2 [{ Sl s 2 } - 1] = -1 + ((si + s 2 - 2) 2 } /2a 2 

[((S1-(S1>) 2 X(S2-(S2>) 2 >J 



< s / a / + "exp(V)^- n ; 2r ; e+ ;, (12) 

The variance cr 2 ^ = ^2 = M2 — 1 °f the GWD is then 



>W 



g+l 
~2b7 



13 3 7 



2g 8g 2 I6g 3 384p 



j + 0{g- b ) 
13) 



for large values of g, as given in Eq. (A8) of Ref. ^3 



Note that b g cancels in the ratio: 



Similarly, the peak of p e (s) (i.e. its mode) occurs at 



(14) 



2b, 



1/2 



e~ e 1 
~ ~ Tg 



(15) 



For the special values of g, the peak positions are tabu- 
lated in Table 1. We further note that s max = 0.96, 0.97, 
and 0.98 occur for g= 6.1, 8.2, and 12.4, respectively. 
When dealing with experimental data, such shifts from 
Smax =1 would be difficult to distinguish. 

While there are some formal justifications for the 
GWD as an optimal description of TWDs for arbitrary A, 
arguably the most convincing argument is a comparison 
of the predicted variance with numerical data generated 
from Monte Carlo simulations. We include in the com- 
parisons some Gaussian approximations (viz. approxima- 
tion schemes which lead to TWDs that are Gaussians), 
alluded to earlier, e.g. in Table 1. For these the dimen- 
sional variance of the TWD must scale like (£) 2 . The 
approximations differ in how the dimensionless part de- 
pends on A. The approximation developed in Grenoble 
by Ihle, Pierre-Louis, and Misbah LJ, [3(| focuses on the 
limit of large g, neglecting the entropic interaction in 
that limit. While the variance a 2 cx A -1 / 2 , the propor- 
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tionality constant is 1.8 times that in the Gruber-Mullins 
case (cf. Table 3). One improve this approximation, es- 
pecially for repulsions that are not extremely strong, by 
including the entropic interaction in an average way. This 
is done by replacing A by 



6 8 1012 15 



icff 



(I)" 



A 



(16) 



(Cf. Eq. (|11|1.) In this modified scheme, a 2 oc g^ 1 . 

The physical meaning of A e g has not been adequately 
discussed heretofore. It corresponds to the full strength 
of the inverse-square repulsion between steps, i.e. the 
modification due to the inclusion of entropic interactions. 
From Eq. ([16]) it is obvious that the contribution of the 
entropic interaction, viz. the difference between the total 
and the energetic interaction, as discussed in conjunction 
with Fig. ([2]), is g/2. This quantity, as noted then, de- 
pends sensitively on A. Note also that, remarkably, the 
ratio of the entropic interaction to the total interaction 
is (g/2) /(g/2) 2 — 2/g; this is the fractional contribution 
that is plotted in Fig. 

Barbier et al. at Saclay [3l|, 0, H(| consider instead 
the m-terrace width relative to m(£), comparing the vari- 
ance for large m to the coefficient of In m expected from 
roughening theory, they conclude that a 2 = 4/(7r 2 g), i.e. 
the same form as the modified Grenoble approximation 
but with a proportionality constant K that is 82% as 
large. This information is summarized in Table 3. Note 
that for the Saclay and Wigner approaches, one must as- 
sume all steps interact with an inverse square repulsion; 
for the others, one can treat either that assumption or 
allow only nearest neighbors to interact energetically. 

In Fig. 2] we compare how well all these theoretical ap- 
proximations account for the behavior extracted from a 
well-controlled numerical experiment, Metropolis Monte 
Carlo simulations [4l| of a TSK model. For this model 
(with A=0), the characteristic distance between close ap- 
proaches of neighboring steps is [42| 

?ycon = W 2 /3/4fc s r=((^)/2) 2 sinh 2 (e/2fc B T), (17) 



which is about 35 at fcsT/e = 1/2 and (£) = 10, used in 
our most extensive calculations. For A > 0, meandering 
is suppressed, making this distance is larger. In most 
cases, we used L y = 2000 ^> y co ii, and iV = 100 steps [43 1. 



We used a standard high-quality random-number gener- 
ator (Ran3 [3]) and averaged over 100 runs using differ- 
ent initial seeds. In these runs the variance reached its 
steady-state value after about 3000 MCS (Monte Carlo 
steps per site); we started "taking data" after 10 4 MCS, 
recording results every 10 MCS until reaching 3xl0 4 MCS 

The excellent agreement between the GWD expression 
and numerical data generated with Monte Carlo simula- 
tions is displayed in Fig. [5] The various predictions of 




FIG. 4: The variance o~ 2 (A) vs. A (and vs. g on the up- 
per abscissa) on a logarithmic scale, plotted for the GWD 
(light solid curve) and for the various Gaussian approxima- 
tions: the modified Grenoble ( short-short-long dashed for 
all steps interacting, short-long dashed for NN step interac- 
tions only), Saclay (short-short-long- long dashed line), and 
Gruber-Mullins (long dashed for all steps interacting, short 
dashed for NN steps only). Monte Carlo data are shown as 
• 's, with statistical errors less than the size of the symbols. 
See text for discussion. [Fig. 2a of Ref. [26j] 



the variance are plotted vs. A. A logarithmic scale is 
used for the horizontal axis so as not to give undue vi- 
sual emphasis to larger values of A nor to blur the region 
of rapid variation for small (but non-vanishing) ^4, for 
which an exact calibration point exists. The Wigner re- 
sult is essentially given by Eq. (TT5)) . Table 3 shows that 
the physical values of A range from near up to the mid 
teens. Some pathology is presumably involved in the rare 
cases in which larger values are observed. There are rel- 
atively few reports of small but non-zero values of A. A 
reason might well be that any of the Gaussian approxi- 
mations manifestly fail in this regime, so that before the 
recognition of the utility of the Wigner distribution, one 
could not deal quantitatively with small A 



45] 



To heighten the contrast, the data in Fig. 0] can be re- 
plotted (in Fig. 2b of Ref. [26|) using as the ordinate 
g(A) x <r 2 (A), so that approximations for whicher 2 cx 

~ — 1/2 

A eS ' become horizontal lines. With such rescaling, 
then, the Saclay and the modified Grenoble (all steps in- 
teracting) predictions appear as lines at ordinates 0.405 
and 0.495, respectively (cf. Table 1). The solid curve rep- 
resenting the GWD rises slowly, from 0.4 at A=l toward 
1/2, capturing a similar rise of the MC data. 

USER'S GUIDE: SYNOPSIS OF FINDINGS 

Refs. [2] and [53] contain several figures showing appli- 
cations, to experimental TWDs, of the perspective dis- 
cussed above. In the following we summarize several spe- 
cific ideas that should be of use in confronting data. 
Interaction strength from variance: From Eqs. |T 
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Vicinal 


T(K) 


o 2 


Q 


A 


W/G 


Pt(110)-(lx2) [45] 


298 




2.2 


0.13 


— 


Cu (19,17,17) [38, 46] 


353 


0.122 


4.1 


2.2 


0.77 


Si(lll) [47J 


1173 


0.11 


3.8 


1.7 


0.96 


Cu(l,l,13) [22, 38] 


348 


0.091 


4.8 


3.0 


1.27 


Cu(ll,7,7) [38, 46] 


306 


0.085 


5.1 


4 


1.37 


Cu(lll) [38, 46J 


313 


0.084 


5.0 


3.6 


1.39 


Cu(lll) [38, 46J 


301 


0.073 


6.0 


6.0 


1.58 


Ag(lOO) 


300 


0.073 


6.4 


6.9 


1.58 


Cu(l,l,19) [22, 38] 


320 


0.070 


6.7 


7.9 


1.64 


Si(lll)-(7x7) [48J 


1100 


0.068 


6.4 


7.0 


1.67 


Si(lll)-(lxl)Br [13] 


853 


0.068 


6.4 


7.0 


1.67 


Si(lll)-Ga [49] 


823 


0.068 


6.6 


7.6 


1.67 


Si(lll)-Al 73 [50] 


1040 


0.058 


7.6 


10.5 


1.85 


Cu(l,l,ll) [31] 


300 


0.053 


8.7 


15 


1.95 


Cu(l,l,13) [22, 38] 


285 


0.044 


10 


20 


2.12 


Pt(lll) [51] 


900 


0.020 


24 


135 


2.59 



TABLE II: Compendium of experiments measuring the vari- 
ances of terrace width distributions of vicinal systems, as of 
the end of the last decade. The estimate of A is obtained from 
the (normalized) variance using Eq. {H}, except for the first- 
row entry, which is based on a direct fit using the 2-parameter 
GWD (cf. 6th item of "User's Guide"). W/G stands for the 
ratio of the estimates of the interaction strength based on 
Wigner and Gruber-Mullins perspectives, Aw/Ag, as given 
in Eq. glj). [Condensed from Ref. [H] 



and (|lip. one can estimate the variance from A, but ex- 
perimentalists usually seek the reverse, measuring the 
variance of the TWD and seeking to extract A. An ex- 
cellent estimate [Hj of the GWD-predicted A from the 
variance, based on series expansion of Eq. (fT5|) , is: 



.4 



1 

Tg 



(a 2 )- 2 - 7(a 2 



27 

T 



35 



(18) 



with all four terms needed to provide a good approxi- 
mation over the full physical range of A. The Gaussian 
methods described earlier essentially use just the first 
term of this expression and adjust the prefactor. 
Gaussian Fits of the GWD: Since TWDs for strong 
repulsions are well described by Gaussians, the GWD is 
well approximated by a Gaussian in this limit. In Ref. 
[3!^ a quantitative assessment is given of how closely the 
two distributions correspond as a function of g. At the 
calibration point (for which an exact solution exists) for 
repulsive interactions (g — 4), the relative difference of 
the standard deviation of a Gaussian fitted to P e (s) from 
the actual standard deviation of this GWD (viz. the 
square root of the second moment of P e (s) about its mean 
of unity) is ~1%, and decreases monotonically with in- 
creasing g. (The x 2 1 however, is poorer with a Gaussian 
than with the GWD.) For this range (g > 4) differences 
between estimates of A obtained from GWD and the var- 
ious Gaussian-fit methods come primarily from different 



philosophies of extracting A from a rather than from dif- 
ferences in the fitting schemes. 

Gaussian approximations assume the peak (mode) is 
at s = 1 (I = (£)); in fact, since peak of the GWD must 
lie below one to achieve a mean of one; the tabulations in 
Table 1 bear this out. Only to the extent that the mode 
is close to unity is a Gaussian approximation reasonable. 
Formulas have been derived [38J indicating the errors in 
fitting g due to errors in the first or zeroth moment of 
the distribution. Another criterion for the validity of the 
Gaussian fitting function is that the TWD not be notice- 
ably asymmetric about its peak, i.e. that it be negligible 
for I > 2(£). 

Analyzing TWD skewness: Unfulfilled hopes: 

When g is too small for the TWD to meet the crite- 
ria for adequate fitting by a Gaussian, we hoped that 
one could obtain reliable estimates of A by analyzing the 
skewness. Although we tried a variety of formulas and 



schemes [27|, |37j, it turned out that in the end a fit to 
the GWD was needed, so that skewness did not offer a 
shortcut to g. 

Correcting estimates based on Gruber-Mullins: 

When dealing with tabulations of data analyzed in the 
traditional way [l[ based on the Gruber-Mullins perspec- 
tive, it is useful to recast Eq. (fTB)) in a form that indicates 
the factor by which the estimate Aw based on GWD ex- 
ceeds the traditional estimate Aq: 



A w /Ag [= A W /A G ] « 3 - 21a 2 



l^ + f.e. (19) 



Since A = Aj3 / '(fcgT) 2 , the ratio of the physical interac- 
tion strengths is the same as that of the dimcnsionlcss 
strengths. 

Alternative: fitting with Gamma distribution: Al- 
though the GWD is a simple, single-parameter function, 
it is not a "canned" distribution and that one needs to in- 
put gamma functions. For the smaller values of g when 
a Gaussian is inappropriate, there is a preprogrammed 
function that is available: the Gamma distribution 



Pr{x) 



T{a)6 a 



(20) 



is serviceable if one recasts the data in terms of 1 as s 2 ; 
one the can identify 9 as b~ x and, more importantly, a 
as (g+l)/2. This approach has not yet been tested with 
actual data. 

Wigner Distribution as a 2-parameter fit: The 

GWDs giving the best fits of experimental TWDs some- 
times have first moments that differ somewhat from the 
first moments of the data, especially in cases termed 
"poor data" [U HH which exhibit a small "hump" at 
large values of s, beyond the mode. Moreover, it may be 
desirable to determine the scaling length (the "effective 
mean," which equals the first moment for ideal GWDs) 
and the variance in a single fitting procedure rather than 



8 



to predetermine this length from the first moment. This 
"refined" scaling implies that the argument of P g should 
be £/£, where I denotes the characteristic length deter- 
mined along with g in a two-parameter least-squares fit 
of the data to a GWD, leading to the replacement [H3 |: 



p e ( s ) -> ((i)/i)p e (s(t)/e) 



((£)/I)P e (£/£). (21) 



In the specific applications to data for copper |28] , (£) /£ 
tends to be greater than unity, typically by several per- 
cent, but it is unclear whether this is true for semicon- 
ductors (or even for other metals). In our Monte Carlo 
simulations [26j |. where we have greater control of pu- 
rity and uniformity than in experiments, the optimal £ is 
essentially identical to (£) : a two-parameter fit is unnec- 
essary. 

Effects of Lattice Discreteness: For large values of 
A, the continuum picture breaks down and one must con- 
front the discreteness of the lattice both in actual physical 
systems and in Monte Carlo simulations. This problem 
is discussed extensively in Refs. [28( and [26j but are not 
repeated here since the relevant values of A are larger 
than found in many physical systems. By adapting the 
celebrated VGL model |53|], we have found [26| that the 
roughening criterion translates to An = (£) 4 /6.Note that 
vicinal surfaces are rough, in contrast to the flat terrace 
orientation. Thus, for (£)=3, e.g., when An exceeds ^14 
(so above the physical range), the vicinal orientation be- 
comes a facet. In some cases such as Si(l 13) , the steps 
are exceedingly straight and uniformly spaced, making 
them excellent templates for growth of nanowires [54j . It 
is therefore likely that this surface is a facet rather than 
a rough vicinal. 

Estimate of Number of Independent Measure- 
ments: In order to estimate uncertainties in the deter- 
mination of the TWD and, ultimately, A, it is important 
to have a realistic value of the number of independent 
measurements, a number generally much smaller than 
the total number of measurements. This problem is dis- 
cussed in Ref. [38[ The upshot is that the number of 
"independent" terrace widths is reduced from L y N (the 
number of atomic spacings across the sampled area along 
the mean step direction times the number of steps) by up 
to nearly two orders of magnitude, due to slow decay of 
correlations perpendicular to the steps. Thus, several 
STM images are needed to obtain decent statistics. 

NEW DEVELOPMENTS AND DIRECTIONS 

Pair correlation function 

The TWD amounts to a many-particle correlation 
function since one insists that there be no step between 
and £ (or s). Instead one can study the step-step correla- 
tion function /i e (s), essentially the probability of finding 
a pair of steps separated by £ regardless of how many 
steps lie in between. Until £ reaches ~ (£) there is little 



difference between the two, but then the pair correla- 
tion begins to rise and peak near 2{£) and at subsequent 
multiples of {£), with a decreasing envelope around the 
oscillations, so that eventually h — » For fitting 

experimental data (in this case, Si(lll) at 1100°C), we 
(Ha | used an asymptotic expression [56| 



(£)h e (s) 
dj(g) 



1 



n 2 g s 2 



£(g) 



= r i 



cos(2ttj's), (22) 



2nm 
Q 



2m 
Q 



that works well for s > 1/2, better than the convention- 
ally used harmonic lattice approximation, since the con- 
ventional "harmonic" approximation [571 ] is insufficiently 
accurate. For smaller s we patched onto a e s e . While 
Ha [Hj| derived a general exact solution for h e (s), it is 
computationally intractable. 
Fokker-Planck 

Recently, Pimpinelli et al. starting from Dyson's 
ID Coulomb gas model and making plausible assump- 
tions, derived a Fokker-Planck equation 



dP(s,t) _ d_ 



2b g s 



§^l P (sM (23) 



where i = t/r, and the characteristic time r is (£) 2 
over the squared strength of the noise in the underlying 
Langevin equation. The steady-state solution of Eq. (|2"5|) 
is the GWD. With it we can describe with compact an- 
alytic expressions the evolution toward equilibrium of 
steps from several experimentally relevant initial distri- 
butions: perfect cleaved crystals (P(s) = S(s — 1) as in 
Fig- Hi step bunches, and prequench equilibrated distri- 
butions at different temperatures (P go (s)). The decay 
time t of the difference of variance of P(s,t) from its 
saturation value, we have also found ^59], can be related 
to underlying atomistic processes in kinetic Monte Carlo 
simulations of the evolution of an initially constrained 
vicinal surface. 

These ideas have broad ramifications. In econophysics 
the same formalism arises for the distribution of the 
means of stock prices in the Heston model [60j]. Similar 
behavior may occur in precipation patterns at geother- 
mal hot springs [611 ]. 

The argument |9j for Eq. ([23|) bolsters the formal justifi- 
cation for the GWD, as does work by Richards et al. [62| 
that uses a two-particle Calogero model [19] (har- moni- 
cally bound interacting spinless fermions on a line). 
Azimuthal Misorientation 

Most of the above has tacitly assumed that the steps 
are close-packed or at least oriented along high symme- 
try directions. For vicinal surfaces surfaces that are mis- 
oriented in the azimuthal in addition to the polar ori- 
entation, there are complications in applying the GWD 
formalism, particularly in determining the dependence of 
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A — ultimately of (3 and A — on in-plane misorientation 9. 
For (3 we have made substantial progress recently, again 
prompted by collaboration with the Jiilich group [63( and 
by persistent probing questions by Ibach. For the {100} 
and {111} faces of fee cubic crystals, we have derived 
surprisingly simple formulas for the (3(9) and (3(9) for 
low T (compared to the terrace roughening temperature) 
64l 65|; in the case of Cu and other noble metals, this cri- 
terion is easily satisfied for room T. For arbitrary T we 
have generated a more complicated analytic expression 
that is straightforward incorporate in continuum simu- 
lations, in particular finite-element codes [66j |. We have 
also carried out ab initio calculations of the characteris- 
tic energies of the lattice-gas models used to understand 
(3(9) for Cu, to insure that these numbers are consistent 
with those estimated from experiments using statistical- 
mechanics reasoning 67J . Fuller discussion is beyond the 
limits of the article. 

Much less is known about A(9), but again the Jiilich 
group is at the forefront of these investigations [HI] . An- 
other important open is the relation between A and sur- 
face stress 69]. There has been no successful extension 
of the seminal theory for isotropic substrates [H, 701 to 
account quantitatively for A of a realistic surface [71 1 . 



CONCLUSION 

The GWD has proved a powerful tool to link the study 
of the step spacings on vicinal surfaces to very general 
ideas of fluctuations. It is remarkable that H. Ibach could 
deduce its form, for the A = case, from physical insight. 
His perspicacious graph (Fig. [3] spawned great progress 
in understanding and exploiting the fluctuations of the 
spacings between steps. Even if one chooses to fit with 
a Gaussian, the theory of GWD provides the most reli- 
able way to extract the strength of the step repulsions 
from the TWD. In addition to clarifying the equilibrium 
properties of steps, these ideas are now helping us to un- 
derstand non-equilibrium aspects, notably the relaxation 
of steps toward equilibrium. Pimpinelli and I are also 
actively investigating applications to various aspects of 
growth in surface and interface problems. 
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Appendix: How the Hamiltonian symmetries lead to 
the Wigner surmise exponents 

This appendix, based on Ref. [72| (see also §4.3 of [25j|). 
describes how the exponents of s in Eqs. ((6H7|) come 
about. For physicists, random matrix theory is rooted 
in the study of energy eigenvalues in nuclear physics. 



As recounted beautifully in Ref. |24|, conventional sta- 
tistical mechanics treats ensembles of identical physical 
systems with the same Hamiltonian but different initial 
conditions. In contrast, Wigner considered ensembles of 
dynamical systems governed by different Hamiltonians 
Ti having some common symmetry property and sought 
generic properties of such ensembles associated with the 
symmetry. Dyson extended this work to show that there 
are three symmetries of the matrices of the Hamilto- 
nian: orthogonal (real symmetric) matrices correspond to 
time-reversal invariance with rotational symmetry, uni- 
tary (complex) matrices to violated time-reversal invari- 
ance (as for electrons in a magnetic field), symplectic 
(see Eq. (2.3) of 0, §2.4 of fl7T|) to time-reversal invari- 
ance with broken rotational symmetry and 1/2-integer 
spin. Wigner then, for convenience, considered Gaussian 
weightings, 



p(H) ~exp[-WVtr(H 2 )], 



(24) 



with the idea that for large matrices the fluctuations of 
the eigenvalues should be independent of the weight fac- 
tors as well as of the specific form of the level spectrum. 
(There is no information about the average values.) 

As described in Ref. [72| ]. we start with the simplest 
case, an orthogonal ensemble (/3=1) with N=2 (just 
2 particles), with TL having diagonal elements h\\ and 
and off-diagonal element /ii2(= ^21 )■ We focus on 
the Jacobean associated with a change of variables for 
the Gaussian-ensemble probability distribution function 
in going from the joint probability distribution function 
p(Ei, E2) for adjacent eigenenergies E\ 1 E 2 to P(s)ds; we 
define variables h = (h\\ + /122) /2, u = h\\ — h 2 2 and s = 
(u 2 +Ah\ 2 ) 1/2 = \E 2 - £i|. Thus, Ei <2 = h±s/2. We 
must now take into account all possible matrix elements 
hn, h%2i h\2- From Eq. (|24[) and with dh\\dh,2i = dhdu, 



p dh 11 dh 22 dh 12 = Jds jj cxp [-2&(£? +£$)] dhdu 



dhyi 



ds 
(25) 

Hence we can identify the inner double integral of 
Eq. (JUD as p(s). Since h u = ±(l/2)(s 2 - u 2 ) 1/2 , 
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\dh 12 /ds\ = (s/2) (s 2 -u 2 ) 1/2 . Since E h2 do not de- 
pend on u, we can pull the exponential out of the inte- 
gration over u, leaving the elementary integral 



(s/2) 



(s 2 -u 2 ) 



-1/2 



du = tts/2 



(26) 



Next, since 2(E\ + E 2 ) = s 2 + 4h 2 , the integration 
over h is also elementary, and we are left with P(s) ~ 
s exp(— bs 2 ). This exact result for N — 2 provides an ex- 
cellent approximation for large N as well. Indeed, near 
the level crossing (corresponding to s — ► 0) the problem 
tends to reduce to a 2 x 2 problem. The integral for P(s) 
over variables u and h\ 2 includes in its integrand a Dirac 
delta function 8 (s — [u 2 + 4ft-f 2 ] 1 ^ 2 )i which vanishes for 
s=0 only when the two [squared] independent variables 
do. Hence, P(s) oc s, corresponding to a circular shell in 
parameter space, and we again have f3=l. 

For unitary ensembles there is an additional indepen- 
dent parameter since h\ 2 becomes (3?e h\ 2 ) 2 + (3m/ii 2 ) 2 . 
Hence, P(s) oc s 2 , corresponding to a spherical shell in 
parameter space, i.e. (3=2. The argument for the sym- 
plectic ensemble leads similarly — but via quaternions or 
Pauli spin matrices — to (3=4. 
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